About
John P drew our attention to various species that were not shown as
common in the NVC standards, but which we frequently find in our sites,
notwithstanding they have been assessed as belonging to a community that
“should” not have them at high frequency. Prominent among these is
Stellaria graminea. The purpose of this notebook is to track
the development of an application to find these unconforming
species.
Start by loading the data
library(plotly)
source("db_extract.R")
the_data <- GetTheData()
And we’re going to need species_frequency by assembly; but maybe
could decide to go by community?
frequency_by_assembly <- FrequencyByAssembly(the_data)
`summarise()` has grouped output by 'assembly_id'. You can override using the `.groups` argument.
`summarise()` has grouped output by 'assembly_id'. You can override using the `.groups` argument.
frequency_by_community <- FrequencyByCommunity(the_data)
`summarise()` has grouped output by 'community', 'species_name'. You can override using the
`.groups` argument.
Let’s plot mean frequency by assembly vs frequency by community for
Stellaria graminea:
fba <- frequency_by_assembly %>% select(assembly_id, species_name, community, freq)
fbc <- frequency_by_community %>% select(species_name,community, freq)
jf <- left_join(fbc, fba, by = c("community", "species_name")) %>% filter(species_name == "Stellaria_graminea" & grepl("MG", community))
f <- ggplot(jf, aes(freq.x, freq.y, colour = community)) +
geom_point()+
xlab("frequency by community") + ylab("frequency by assembly")+
geom_segment(aes(x = 0, xend = 1, y = 0, yend = 1, colour ="red"))
print(f)

This is for S. graminea and the MG communities. So (a) there
are lots more assemblies than there are communities (of course) and (b)
frequency by assembly tends to be larger than frequency by community;
aggregating over community has smoothed out a lot of variation. My
feeling is that it will be more illuminating to eliminate the assembly
variable to start with, we can explore variability with assembly
later.
Next question: how to compare with standard values. Explore select
into table mg_rodwell in the meadows database … in which we have columns
community, species_id (and name, but not to be relied on) and p_central,
the central frequency of the categories I .. V.
Stellaria graminea id 139. It so happens that p_central ==
0.1 for S. graminea for all MG communities: SELECT Community,
species_id, p_central FROM meadows.mg_rodwell where Community like “MG%”
&& species_id = 139; ( I wonder whether I can communicate
directly with the DB from here using SQL?).
Anyway, it seems quite possible to collect the p_central values for
all species and plot against frequency by community
q <- "SELECT Community, species_id, p_central FROM meadows.mg_rodwell where Community like 'MG%';"
std_freqs <- query(q)
Warning in .local(conn, statement, ...) :
Decimal MySQL column 2 imported as numeric
and match up with species frequencies by community
rm(jf) # We don't need it any more
jf1 <- left_join(frequency_by_community, std_freqs, by = c("community"="Community", "species_id"))
Lots of NAs here; (a) because mires (M) are included in the survey
data but not in the mg_rodwell table - yet - because as its name
implies, it is just the MG communities. And (b), because it seems we
have detected some species that the standards don’t list in those
communities at all. So (a), filter out all the mires; (b) set the
remaining NAs to zero.
# Would be easier if we'd excluded mires from the_data to start with!
jf2 <- jf1 %>% filter(grepl("MG", community)) %>% replace_na(list(p_central = 0))
rm(jf1)
So now we want filter to include only cases where CrI5 > p_central
or CrI95 < p_central, and then plot (survey) frequencies against
standard frequencies (p_central)
jf3 <- jf2 %>% filter(CrI5 > p_central || CrI95 < p_central)
f2 <- ggplot(jf3, aes(p_central, freq, colour = community)) +
geom_jitter(aes(text = species_name), size = 3) + #, name = "species_name")+
xlab("Standard frequency") + ylab("Survey frequency")+
geom_segment(aes(x = 0, xend = 1, y = 0, yend = 1, colour ="red"))
Warning: Ignoring unknown aesthetics: text
ggplotly(f2) %>% layout(legend = list(orientation = "h", x = 0.4, y = -0.2))
NA
This gives a surprising number of species in strong disagrrement with
the standard (we may ignore the big cluster near (0, 0) because of the
way CrI5 is calculated - it can be very near zero, but not actually
zero). Interactive plot: is this useful?
Even so, there’s too many species here for this to be a useful way of
exploring the data. Filter out low frequency species, and then perhaps
present the data as an interactive table.
Reduce species count.
There are 22991 records. Stellaria graminea accounts for 540
or 2.3487452%, so suggest reducing the data to species with less than 2%
of the count. To start with, we have 300 species.
d1 <- (the_data
%>% select(records_id, species_id)
%>% group_by(species_id)
%>% summarise(cnt = n())
%>% mutate(frac = cnt/nrow(the_data))
%>% filter(frac >= 0.02))
Removing those accounting for less than 2% of the records leaves 13
species. Make the reduced dataset and plot survey frequency vs standard
again, coded by species:
reduced <- left_join(d1, the_data, by = "species_id")
jf4 <- left_join(d1, jf3, by = "species_id")
f3 <- ggplot(jf4, aes(p_central, freq, colour = species_name)) +
geom_jitter(aes(text = community), size = 3) +
xlab("Standard frequency") + ylab("Survey frequency")+
geom_segment(aes(x = 0, xend = 1, y = 0, yend = 1), colour ="red")
Warning: Ignoring unknown aesthetics: text
ggplotly(f3) %>% layout(legend = list(orientation = "h", x = 0.4, y = -0.2))
NA
This looks much more useful. Community is now available in the
interactive part. Using 13 classes: on reflection (a) it might be better
to use the top 12 classes (allowing to use ColorBrewer Paired Class 12
for the legend); and (b), ultimately to analyse by community (group, as
used for poster BES2019) and select the top twelve species for each
community.
Doing (b) first: sort out mg5a data:
cf <- head(the_data
%>% filter(community == "MG5a")
%>% group_by(species_id)
%>% summarise(hits = n()) %>% arrange(desc(hits)), 12)
d_mg5a <- the_data %>% filter(community == "MG5a") %>% right_join(cf, by = "species_id")
fc_mg5a <- FrequencyByCommunity(d_mg5a)
`summarise()` has grouped output by 'community', 'species_name'. You can override using the
`.groups` argument.
# Add standard frequencies
fc_mg5a <- (left_join(fc_mg5a, std_freqs, by = c("community"="Community", "species_id"))
%>% mutate(std_freq = replace_na(p_central, 0))
%>% select(-p_central))
f4 <- ggplot(fc_mg5a, aes(std_freq, freq, colour = species_name)) +
geom_pointrange(aes(ymin = CrI5, ymax = CrI95, text = community), size = 3) +
scale_colour_brewer(palette="Paired") +
xlab("Standard frequency") + ylab("Survey frequency") +
geom_segment(aes(x = 0, xend = 1, y = 0, yend = 1), colour ="red")
Warning: Ignoring unknown aesthetics: text
ggplotly(f4) %>% layout(legend = list(orientation = "h", x = 0.4, y = -0.2))
NA
The credibility intervals are tiny because we are working with 200 -
300 hits and over 400 trials. It will be interesting to see how this
changes on an assembly basis, and in communities where we have less
data.
Let’s try the same thing but on an assembly basis:
fa_mg5a <- FrequencyByAssembly(d_mg5a)
`summarise()` has grouped output by 'assembly_id'. You can override using the `.groups` argument.
`summarise()` has grouped output by 'assembly_id'. You can override using the `.groups` argument.
# Add standard frequencies
fa_mg5a <- (left_join(fa_mg5a, std_freqs, by = c("community"="Community", "species_id"))
%>% filter(!is.na(species_id))
%>% filter(!is.na(community))
%>% mutate(std_freq = replace_na(p_central, 0))
%>% select(-p_central))
pd <- position_dodge(0.2)
f5 <- ggplot(fa_mg5a, aes(std_freq, freq, colour = species_name)) +
geom_pointrange(aes(text = paste(assembly_name, community), ymin = CrI5, ymax = CrI95), size = 1, position = pd) +
scale_colour_brewer(palette="Paired") +
xlab("Standard frequency") + ylab("Survey frequency") +
geom_segment(aes(x = 0, xend = 1, y = 0, yend = 1), colour ="red")
Warning: Ignoring unknown aesthetics: text
ggplotly(f5) %>% layout(legend = list(orientation = "h", x = 0.4, y = -0.2))
NA
Which certainly makes the point about the credibility intervals, but
is not otherwise very useful, so have to think of another way to deal
with assembly-level analysis.
Single assembly is nice:
f6 <- ggplot(head(fa_mg5a, 11), aes(std_freq, freq, colour = species_name)) +
geom_pointrange(aes(text = assembly_name, ymin = CrI5, ymax = CrI95), size = 1, position = pd) +
scale_colour_brewer(palette="Paired") +
xlab("Standard frequency") + ylab("Survey frequency") +
geom_segment(aes(x = 0, xend = 1, y = 0, yend = 1), colour ="red")
Warning: Ignoring unknown aesthetics: text
ggplotly(f6) %>% layout(legend = list(orientation = "h", x = 0.4, y = -0.2))
NA
but really needs interactivity so we can search assembly-by-assembly.
Bear in mind we have total 251 assemblies, of which 75 are MG5a.
Looks like this will need special treatment with a Shiny app, beyond
the scope of this exploratory notebook.
LS0tDQp0aXRsZTogIlN0ZWxsYXJpYSBub3RlYm9vayINCm91dHB1dDogaHRtbF9ub3RlYm9vaw0KZmlnLndpZHRoOiA2DQpmaWcuaGVpZ2h0OiA0DQotLS0NCiMgQWJvdXQNCg0KSm9obiBQIGRyZXcgb3VyIGF0dGVudGlvbiB0byB2YXJpb3VzIHNwZWNpZXMgdGhhdCB3ZXJlIG5vdCBzaG93biBhcyBjb21tb24gaW4gdGhlIE5WQyBzdGFuZGFyZHMsIGJ1dCB3aGljaCB3ZSBmcmVxdWVudGx5IGZpbmQgaW4gb3VyIHNpdGVzLCBub3R3aXRoc3RhbmRpbmcgdGhleSBoYXZlIGJlZW4gYXNzZXNzZWQgYXMgYmVsb25naW5nIHRvIGEgY29tbXVuaXR5IHRoYXQgInNob3VsZCIgbm90IGhhdmUgdGhlbSBhdCBoaWdoIGZyZXF1ZW5jeS4gUHJvbWluZW50IGFtb25nIHRoZXNlIGlzICpTdGVsbGFyaWEgZ3JhbWluZWEqLiBUaGUgcHVycG9zZSBvZiB0aGlzIG5vdGVib29rIGlzIHRvIHRyYWNrIHRoZSBkZXZlbG9wbWVudCBvZiBhbiBhcHBsaWNhdGlvbiB0byBmaW5kIHRoZXNlIHVuY29uZm9ybWluZyBzcGVjaWVzLg0KDQpTdGFydCBieSBsb2FkaW5nIHRoZSBkYXRhDQoNCmBgYHtyfQ0KbGlicmFyeShwbG90bHkpDQpzb3VyY2UoImRiX2V4dHJhY3QuUiIpDQp0aGVfZGF0YSA8LSBHZXRUaGVEYXRhKCkNCmBgYA0KDQpBbmQgd2UncmUgZ29pbmcgdG8gbmVlZCBzcGVjaWVzX2ZyZXF1ZW5jeSBieSBhc3NlbWJseTsgYnV0IG1heWJlIGNvdWxkIGRlY2lkZSB0byBnbyBieSBjb21tdW5pdHk/DQpgYGB7cn0NCmZyZXF1ZW5jeV9ieV9hc3NlbWJseSA8LSBGcmVxdWVuY3lCeUFzc2VtYmx5KHRoZV9kYXRhKQ0KZnJlcXVlbmN5X2J5X2NvbW11bml0eSA8LSBGcmVxdWVuY3lCeUNvbW11bml0eSh0aGVfZGF0YSkNCmBgYA0KTGV0J3MgcGxvdCBtZWFuIGZyZXF1ZW5jeSBieSBhc3NlbWJseSB2cyBmcmVxdWVuY3kgYnkgY29tbXVuaXR5IGZvciAqU3RlbGxhcmlhIGdyYW1pbmVhKjoNCmBgYHtyfQ0KZmJhIDwtIGZyZXF1ZW5jeV9ieV9hc3NlbWJseSAlPiUgc2VsZWN0KGFzc2VtYmx5X2lkLCBzcGVjaWVzX25hbWUsIGNvbW11bml0eSwgZnJlcSkNCmZiYyA8LSBmcmVxdWVuY3lfYnlfY29tbXVuaXR5ICU+JSBzZWxlY3Qoc3BlY2llc19uYW1lLGNvbW11bml0eSwgZnJlcSkNCmpmIDwtIGxlZnRfam9pbihmYmMsIGZiYSwgYnkgPSBjKCJjb21tdW5pdHkiLCAic3BlY2llc19uYW1lIikpICU+JSBmaWx0ZXIoc3BlY2llc19uYW1lID09ICJTdGVsbGFyaWFfZ3JhbWluZWEiICYgZ3JlcGwoIk1HIiwgY29tbXVuaXR5KSkNCmYgPC0gZ2dwbG90KGpmLCBhZXMoZnJlcS54LCBmcmVxLnksIGNvbG91ciA9IGNvbW11bml0eSkpICsNCiAgZ2VvbV9wb2ludCgpKw0KICB4bGFiKCJmcmVxdWVuY3kgYnkgY29tbXVuaXR5IikgKyB5bGFiKCJmcmVxdWVuY3kgYnkgYXNzZW1ibHkiKSsNCiAgZ2VvbV9zZWdtZW50KGFlcyh4ID0gMCwgeGVuZCA9IDEsIHkgPSAwLCB5ZW5kID0gMSwgY29sb3VyID0icmVkIikpDQpwcmludChmKQ0KYGBgDQoNCg0KYGBge3J9DQpgYGANClRoaXMgaXMgZm9yICpTLiBncmFtaW5lYSogYW5kIHRoZSBNRyBjb21tdW5pdGllcy4gU28gKGEpIHRoZXJlIGFyZSBsb3RzIG1vcmUgYXNzZW1ibGllcyB0aGFuIHRoZXJlIGFyZSBjb21tdW5pdGllcyAob2YgY291cnNlKSBhbmQgKGIpIGZyZXF1ZW5jeSBieSBhc3NlbWJseSB0ZW5kcyB0byBiZSBsYXJnZXIgdGhhbiBmcmVxdWVuY3kgYnkgY29tbXVuaXR5OyBhZ2dyZWdhdGluZyBvdmVyIGNvbW11bml0eSBoYXMgc21vb3RoZWQgb3V0IGEgbG90IG9mIHZhcmlhdGlvbi4gTXkgZmVlbGluZyBpcyB0aGF0IGl0IHdpbGwgYmUgbW9yZSBpbGx1bWluYXRpbmcgdG8gZWxpbWluYXRlIHRoZSBhc3NlbWJseSB2YXJpYWJsZSB0byBzdGFydCB3aXRoLCB3ZSBjYW4gZXhwbG9yZSB2YXJpYWJpbGl0eSB3aXRoIGFzc2VtYmx5IGxhdGVyLg0KDQpOZXh0IHF1ZXN0aW9uOiBob3cgdG8gY29tcGFyZSB3aXRoIHN0YW5kYXJkIHZhbHVlcy4gRXhwbG9yZSBzZWxlY3QgaW50byB0YWJsZSBtZ19yb2R3ZWxsIGluIHRoZSBtZWFkb3dzIGRhdGFiYXNlIC4uLiBpbiB3aGljaCB3ZSBoYXZlIGNvbHVtbnMgY29tbXVuaXR5LCBzcGVjaWVzX2lkIChhbmQgbmFtZSwgYnV0IG5vdCB0byBiZSByZWxpZWQgb24pIGFuZCBwX2NlbnRyYWwsIHRoZSBjZW50cmFsIGZyZXF1ZW5jeSBvZiB0aGUgY2F0ZWdvcmllcyBJIC4uIFYuDQoNCipTdGVsbGFyaWEgZ3JhbWluZWEqIGlkIDEzOS4gSXQgc28gaGFwcGVucyB0aGF0IHBfY2VudHJhbCA9PSAwLjEgZm9yICpTLiBncmFtaW5lYSogZm9yIGFsbCBNRyBjb21tdW5pdGllczoNClNFTEVDVCBDb21tdW5pdHksIHNwZWNpZXNfaWQsIHBfY2VudHJhbCBGUk9NIG1lYWRvd3MubWdfcm9kd2VsbCB3aGVyZSBDb21tdW5pdHkgbGlrZSAiTUclIiAmJiBzcGVjaWVzX2lkID0gMTM5Ow0KKCBJIHdvbmRlciB3aGV0aGVyIEkgY2FuIGNvbW11bmljYXRlIGRpcmVjdGx5IHdpdGggdGhlIERCIGZyb20gaGVyZSB1c2luZyBTUUw/KS4NCg0KQW55d2F5LCBpdCBzZWVtcyBxdWl0ZSBwb3NzaWJsZSB0byBjb2xsZWN0IHRoZSBwX2NlbnRyYWwgdmFsdWVzIGZvciBhbGwgc3BlY2llcyBhbmQgcGxvdCBhZ2FpbnN0IGZyZXF1ZW5jeSBieSBjb21tdW5pdHkNCg0KYGBge3J9DQpxIDwtICJTRUxFQ1QgQ29tbXVuaXR5LCBzcGVjaWVzX2lkLCBwX2NlbnRyYWwgRlJPTSBtZWFkb3dzLm1nX3JvZHdlbGwgd2hlcmUgQ29tbXVuaXR5IGxpa2UgJ01HJSc7Ig0Kc3RkX2ZyZXFzIDwtIHF1ZXJ5KHEpDQpgYGANCmFuZCBtYXRjaCB1cCB3aXRoIHNwZWNpZXMgZnJlcXVlbmNpZXMgYnkgY29tbXVuaXR5DQpgYGB7cn0NCnJtKGpmKSAjIFdlIGRvbid0IG5lZWQgaXQgYW55IG1vcmUNCmpmMSA8LSBsZWZ0X2pvaW4oZnJlcXVlbmN5X2J5X2NvbW11bml0eSwgc3RkX2ZyZXFzLCBieSA9IGMoImNvbW11bml0eSI9IkNvbW11bml0eSIsICJzcGVjaWVzX2lkIikpDQoNCmBgYA0KDQpMb3RzIG9mIE5BcyBoZXJlOyAoYSkgYmVjYXVzZSBtaXJlcyAoTSkgYXJlIGluY2x1ZGVkIGluIHRoZSBzdXJ2ZXkgZGF0YSBidXQgbm90IGluIHRoZSBtZ19yb2R3ZWxsIHRhYmxlIC0geWV0IC0gYmVjYXVzZSBhcyBpdHMgbmFtZSBpbXBsaWVzLCBpdCBpcyBqdXN0IHRoZSBNRyBjb21tdW5pdGllcy4gQW5kIChiKSwgYmVjYXVzZSBpdCBzZWVtcyB3ZSBoYXZlIGRldGVjdGVkIHNvbWUgc3BlY2llcyB0aGF0IHRoZSBzdGFuZGFyZHMgZG9uJ3QgbGlzdCBpbiB0aG9zZSBjb21tdW5pdGllcyBhdCBhbGwuIFNvIChhKSwgZmlsdGVyIG91dCBhbGwgdGhlIG1pcmVzOyAoYikgc2V0IHRoZSByZW1haW5pbmcgTkFzIHRvIHplcm8uDQpgYGB7cn0NCiMgV291bGQgYmUgZWFzaWVyIGlmIHdlJ2QgZXhjbHVkZWQgbWlyZXMgZnJvbSB0aGVfZGF0YSB0byBzdGFydCB3aXRoISANCmpmMiA8LSBqZjEgJT4lIGZpbHRlcihncmVwbCgiTUciLCBjb21tdW5pdHkpKSAlPiUgcmVwbGFjZV9uYShsaXN0KHBfY2VudHJhbCA9IDApKQ0Kcm0oamYxKQ0KYGBgDQpTbyBub3cgd2Ugd2FudCBmaWx0ZXIgdG8gaW5jbHVkZSBvbmx5IGNhc2VzIHdoZXJlIENySTUgPiBwX2NlbnRyYWwgb3IgQ3JJOTUgPCBwX2NlbnRyYWwsIGFuZCB0aGVuIHBsb3QgKHN1cnZleSkgZnJlcXVlbmNpZXMgYWdhaW5zdCBzdGFuZGFyZCBmcmVxdWVuY2llcyAocF9jZW50cmFsKQ0KYGBge3J9DQpqZjMgPC0gamYyICU+JSBmaWx0ZXIoQ3JJNSA+IHBfY2VudHJhbCB8fCBDckk5NSA8IHBfY2VudHJhbCkNCmYyIDwtIGdncGxvdChqZjMsIGFlcyhwX2NlbnRyYWwsIGZyZXEsIGNvbG91ciA9IGNvbW11bml0eSkpICsNCiAgZ2VvbV9qaXR0ZXIoYWVzKHRleHQgPSBzcGVjaWVzX25hbWUpLCBzaXplID0gMykgKyAjLCBuYW1lID0gInNwZWNpZXNfbmFtZSIpKw0KICB4bGFiKCJTdGFuZGFyZCBmcmVxdWVuY3kiKSArIHlsYWIoIlN1cnZleSBmcmVxdWVuY3kiKSsNCiAgZ2VvbV9zZWdtZW50KGFlcyh4ID0gMCwgeGVuZCA9IDEsIHkgPSAwLCB5ZW5kID0gMSwgY29sb3VyID0icmVkIikpDQpnZ3Bsb3RseShmMikgJT4lIGxheW91dChsZWdlbmQgPSBsaXN0KG9yaWVudGF0aW9uID0gImgiLCB4ID0gMC40LCB5ID0gLTAuMikpDQoNCmBgYA0KVGhpcyBnaXZlcyBhIHN1cnByaXNpbmcgbnVtYmVyIG9mIHNwZWNpZXMgaW4gc3Ryb25nIGRpc2FncnJlbWVudCB3aXRoIHRoZSBzdGFuZGFyZCAod2UgbWF5IGlnbm9yZSB0aGUgYmlnIGNsdXN0ZXIgbmVhciAoMCwgMCkgYmVjYXVzZSBvZiB0aGUgd2F5IENySTUgaXMgY2FsY3VsYXRlZCAtIGl0IGNhbiBiZSB2ZXJ5IG5lYXIgemVybywgYnV0IG5vdCBhY3R1YWxseSB6ZXJvKS4gSW50ZXJhY3RpdmUgcGxvdDogaXMgdGhpcyB1c2VmdWw/DQoNCkV2ZW4gc28sIHRoZXJlJ3MgdG9vIG1hbnkgc3BlY2llcyBoZXJlIGZvciB0aGlzIHRvIGJlIGEgdXNlZnVsIHdheSBvZiBleHBsb3JpbmcgdGhlIGRhdGEuIEZpbHRlciBvdXQgbG93IGZyZXF1ZW5jeSBzcGVjaWVzLCBhbmQgdGhlbiBwZXJoYXBzIHByZXNlbnQgdGhlIGRhdGEgYXMgYW4gaW50ZXJhY3RpdmUgdGFibGUuDQoNCiMjIyBSZWR1Y2Ugc3BlY2llcyBjb3VudC4NClRoZXJlIGFyZSBgciBucm93KHRoZV9kYXRhKWAgcmVjb3Jkcy4gKlN0ZWxsYXJpYSBncmFtaW5lYSogYWNjb3VudHMgZm9yIGByIHRoZV9kYXRhICU+JSBzZWxlY3Qoc3BlY2llc19pZCwgcmVjb3Jkc19pZCkgJT4lIGZpbHRlcihzcGVjaWVzX2lkID09IDEzOSkgJT4lIG5yb3coKWANCm9yIGByIDEwMCoodGhlX2RhdGEgJT4lIHNlbGVjdChzcGVjaWVzX2lkLCByZWNvcmRzX2lkKSAlPiUgZmlsdGVyKHNwZWNpZXNfaWQgPT0gMTM5KSAlPiUgbnJvdygpKS9ucm93KHRoZV9kYXRhKWAlLCBzbyBzdWdnZXN0IHJlZHVjaW5nIHRoZSBkYXRhIHRvIHNwZWNpZXMgd2l0aCBsZXNzIHRoYW4gMiUgb2YgdGhlIGNvdW50LiBUbyBzdGFydCB3aXRoLCB3ZSBoYXZlIGByIHRoZV9kYXRhICU+JSBzZWxlY3Qoc3BlY2llc19pZCkgJT4lIGRpc3RpbmN0KCkgJT4lIG5yb3coKWAgc3BlY2llcy4gDQpgYGB7cn0NCmQxIDwtICh0aGVfZGF0YSANCiAgICAgICAlPiUgc2VsZWN0KHJlY29yZHNfaWQsIHNwZWNpZXNfaWQpIA0KICAgICAgICU+JSBncm91cF9ieShzcGVjaWVzX2lkKSANCiAgICAgICAlPiUgc3VtbWFyaXNlKGNudCA9IG4oKSkNCiAgICAgICAlPiUgbXV0YXRlKGZyYWMgPSBjbnQvbnJvdyh0aGVfZGF0YSkpDQogICAgICAgJT4lIGZpbHRlcihmcmFjID49IDAuMDIpKQ0KYGBgDQpSZW1vdmluZyB0aG9zZSBhY2NvdW50aW5nIGZvciBsZXNzIHRoYW4gMiUgb2YgdGhlIHJlY29yZHMgbGVhdmVzIGByIG5yb3coZDEpYCBzcGVjaWVzLg0KTWFrZSB0aGUgcmVkdWNlZCBkYXRhc2V0IGFuZCBwbG90IHN1cnZleSBmcmVxdWVuY3kgdnMgc3RhbmRhcmQgYWdhaW4sIGNvZGVkIGJ5IHNwZWNpZXM6DQpgYGB7cn0NCnJlZHVjZWQgPC0gbGVmdF9qb2luKGQxLCB0aGVfZGF0YSwgYnkgPSAic3BlY2llc19pZCIpDQpqZjQgPC0gbGVmdF9qb2luKGQxLCBqZjMsIGJ5ID0gInNwZWNpZXNfaWQiKQ0KZjMgPC0gZ2dwbG90KGpmNCwgYWVzKHBfY2VudHJhbCwgZnJlcSwgY29sb3VyID0gc3BlY2llc19uYW1lKSkgKw0KICBnZW9tX2ppdHRlcihhZXModGV4dCA9IGNvbW11bml0eSksIHNpemUgPSAzKSArIA0KICB4bGFiKCJTdGFuZGFyZCBmcmVxdWVuY3kiKSArIHlsYWIoIlN1cnZleSBmcmVxdWVuY3kiKSsNCiAgZ2VvbV9zZWdtZW50KGFlcyh4ID0gMCwgeGVuZCA9IDEsIHkgPSAwLCB5ZW5kID0gMSksIGNvbG91ciA9InJlZCIpDQpnZ3Bsb3RseShmMykgJT4lIGxheW91dChsZWdlbmQgPSBsaXN0KG9yaWVudGF0aW9uID0gImgiLCB4ID0gMC40LCB5ID0gLTAuMikpDQoNCmBgYA0KVGhpcyBsb29rcyBtdWNoIG1vcmUgdXNlZnVsLiBDb21tdW5pdHkgaXMgbm93IGF2YWlsYWJsZSBpbiB0aGUgaW50ZXJhY3RpdmUgcGFydC4gVXNpbmcgMTMgY2xhc3Nlczogb24gcmVmbGVjdGlvbiAoYSkgaXQgbWlnaHQgYmUgYmV0dGVyIHRvIHVzZSB0aGUgdG9wIDEyIGNsYXNzZXMgKGFsbG93aW5nIHRvIHVzZSBDb2xvckJyZXdlciBQYWlyZWQgQ2xhc3MgMTIgZm9yIHRoZSBsZWdlbmQpOyBhbmQgKGIpLCB1bHRpbWF0ZWx5IHRvIGFuYWx5c2UgYnkgY29tbXVuaXR5IChncm91cCwgYXMgdXNlZCBmb3IgcG9zdGVyIEJFUzIwMTkpIGFuZCBzZWxlY3QgdGhlIHRvcCB0d2VsdmUgc3BlY2llcyBmb3IgZWFjaCBjb21tdW5pdHkuDQoNCkRvaW5nIChiKSBmaXJzdDogc29ydCBvdXQgbWc1YSBkYXRhOg0KYGBge3J9DQoNCmNmIDwtIGhlYWQodGhlX2RhdGEgDQogICAgICAgICAgICU+JSBmaWx0ZXIoY29tbXVuaXR5ID09ICJNRzVhIikNCiAgICAgICAgICAgJT4lIGdyb3VwX2J5KHNwZWNpZXNfaWQpIA0KICAgICAgICAgICAlPiUgc3VtbWFyaXNlKGhpdHMgPSBuKCkpICU+JSBhcnJhbmdlKGRlc2MoaGl0cykpLCAxMikNCmRfbWc1YSA8LSB0aGVfZGF0YSAlPiUgZmlsdGVyKGNvbW11bml0eSA9PSAiTUc1YSIpICU+JSByaWdodF9qb2luKGNmLCBieSA9ICJzcGVjaWVzX2lkIikNCmZjX21nNWEgPC0gRnJlcXVlbmN5QnlDb21tdW5pdHkoZF9tZzVhKQ0KIyBBZGQgc3RhbmRhcmQgZnJlcXVlbmNpZXMNCmZjX21nNWEgPC0gKGxlZnRfam9pbihmY19tZzVhLCBzdGRfZnJlcXMsIGJ5ID0gYygiY29tbXVuaXR5Ij0iQ29tbXVuaXR5IiwgInNwZWNpZXNfaWQiKSkgDQogICAgICAgICAgICAlPiUgbXV0YXRlKHN0ZF9mcmVxID0gcmVwbGFjZV9uYShwX2NlbnRyYWwsIDApKSANCiAgICAgICAgICAgICU+JSBzZWxlY3QoLXBfY2VudHJhbCkpDQoNCmY0IDwtIGdncGxvdChmY19tZzVhLCBhZXMoc3RkX2ZyZXEsIGZyZXEsIGNvbG91ciA9IHNwZWNpZXNfbmFtZSkpICsNCiAgZ2VvbV9wb2ludHJhbmdlKGFlcyh5bWluID0gQ3JJNSwgeW1heCA9IENySTk1LCB0ZXh0ID0gY29tbXVuaXR5KSwgc2l6ZSA9IDMpICsgDQogIHNjYWxlX2NvbG91cl9icmV3ZXIocGFsZXR0ZT0iUGFpcmVkIikgKw0KICB4bGFiKCJTdGFuZGFyZCBmcmVxdWVuY3kiKSArIHlsYWIoIlN1cnZleSBmcmVxdWVuY3kiKSArDQogIGdlb21fc2VnbWVudChhZXMoeCA9IDAsIHhlbmQgPSAxLCB5ID0gMCwgeWVuZCA9IDEpLCBjb2xvdXIgPSJyZWQiKQ0KZ2dwbG90bHkoZjQpICU+JSBsYXlvdXQobGVnZW5kID0gbGlzdChvcmllbnRhdGlvbiA9ICJoIiwgeCA9IDAuNCwgeSA9IC0wLjIpKQ0KDQpgYGANCg0KVGhlIGNyZWRpYmlsaXR5IGludGVydmFscyBhcmUgdGlueSBiZWNhdXNlIHdlIGFyZSB3b3JraW5nIHdpdGggMjAwIC0gMzAwIGhpdHMgYW5kIG92ZXIgNDAwIHRyaWFscy4gSXQgd2lsbCBiZSBpbnRlcmVzdGluZyB0byBzZWUgaG93IHRoaXMgY2hhbmdlcyBvbiBhbiBhc3NlbWJseSBiYXNpcywgYW5kIGluIGNvbW11bml0aWVzIHdoZXJlIHdlIGhhdmUgbGVzcyBkYXRhLg0KDQpMZXQncyB0cnkgdGhlIHNhbWUgdGhpbmcgYnV0IG9uIGFuIGFzc2VtYmx5IGJhc2lzOg0KYGBge3J9DQpmYV9tZzVhIDwtIEZyZXF1ZW5jeUJ5QXNzZW1ibHkoZF9tZzVhKQ0KIyBBZGQgc3RhbmRhcmQgZnJlcXVlbmNpZXMNCmZhX21nNWEgPC0gKGxlZnRfam9pbihmYV9tZzVhLCBzdGRfZnJlcXMsIGJ5ID0gYygiY29tbXVuaXR5Ij0iQ29tbXVuaXR5IiwgInNwZWNpZXNfaWQiKSkNCiAgICAgICAgICAgICU+JSBmaWx0ZXIoIWlzLm5hKHNwZWNpZXNfaWQpKQ0KICAgICAgICAgICAgJT4lIGZpbHRlcighaXMubmEoY29tbXVuaXR5KSkNCiAgICAgICAgICAgICU+JSBtdXRhdGUoc3RkX2ZyZXEgPSByZXBsYWNlX25hKHBfY2VudHJhbCwgMCkpIA0KICAgICAgICAgICAgJT4lIHNlbGVjdCgtcF9jZW50cmFsKSkNCg0KcGQgPC0gcG9zaXRpb25fZG9kZ2UoMC4yKQ0KZjUgPC0gZ2dwbG90KGZhX21nNWEsIGFlcyhzdGRfZnJlcSwgZnJlcSwgY29sb3VyID0gc3BlY2llc19uYW1lKSkgKw0KICBnZW9tX3BvaW50cmFuZ2UoYWVzKHRleHQgPSBwYXN0ZShhc3NlbWJseV9uYW1lLCBjb21tdW5pdHkpLCB5bWluID0gQ3JJNSwgeW1heCA9IENySTk1KSwgc2l6ZSA9IDEsIHBvc2l0aW9uID0gcGQpICsNCiAgc2NhbGVfY29sb3VyX2JyZXdlcihwYWxldHRlPSJQYWlyZWQiKSArDQogIHhsYWIoIlN0YW5kYXJkIGZyZXF1ZW5jeSIpICsgeWxhYigiU3VydmV5IGZyZXF1ZW5jeSIpICsNCiAgZ2VvbV9zZWdtZW50KGFlcyh4ID0gMCwgeGVuZCA9IDEsIHkgPSAwLCB5ZW5kID0gMSksIGNvbG91ciA9InJlZCIpICANCmdncGxvdGx5KGY1KSAlPiUgbGF5b3V0KGxlZ2VuZCA9IGxpc3Qob3JpZW50YXRpb24gPSAiaCIsIHggPSAwLjQsIHkgPSAtMC4yKSkNCg0KYGBgDQpXaGljaCBjZXJ0YWlubHkgbWFrZXMgdGhlIHBvaW50IGFib3V0IHRoZSBjcmVkaWJpbGl0eSBpbnRlcnZhbHMsIGJ1dCBpcyBub3Qgb3RoZXJ3aXNlIHZlcnkgdXNlZnVsLCBzbyBoYXZlIHRvIHRoaW5rIG9mIGFub3RoZXIgd2F5IHRvIGRlYWwgd2l0aCBhc3NlbWJseS1sZXZlbCBhbmFseXNpcy4NCg0KU2luZ2xlIGFzc2VtYmx5IGlzIG5pY2U6DQpgYGB7cn0NCmY2IDwtIGdncGxvdChoZWFkKGZhX21nNWEsIDExKSwgYWVzKHN0ZF9mcmVxLCBmcmVxLCBjb2xvdXIgPSBzcGVjaWVzX25hbWUpKSArDQogIGdlb21fcG9pbnRyYW5nZShhZXModGV4dCA9IGFzc2VtYmx5X25hbWUsIHltaW4gPSBDckk1LCB5bWF4ID0gQ3JJOTUpLCBzaXplID0gMSwgcG9zaXRpb24gPSBwZCkgKw0KICBzY2FsZV9jb2xvdXJfYnJld2VyKHBhbGV0dGU9IlBhaXJlZCIpICsNCiAgeGxhYigiU3RhbmRhcmQgZnJlcXVlbmN5IikgKyB5bGFiKCJTdXJ2ZXkgZnJlcXVlbmN5IikgKw0KICBnZW9tX3NlZ21lbnQoYWVzKHggPSAwLCB4ZW5kID0gMSwgeSA9IDAsIHllbmQgPSAxKSwgY29sb3VyID0icmVkIikNCiAgZ2dwbG90bHkoZjYpICU+JSBsYXlvdXQobGVnZW5kID0gbGlzdChvcmllbnRhdGlvbiA9ICJoIiwgeCA9IDAuNCwgeSA9IC0wLjIpKQ0KDQpgYGANCg0KYnV0IHJlYWxseSBuZWVkcyBpbnRlcmFjdGl2aXR5IHNvIHdlIGNhbiBzZWFyY2ggYXNzZW1ibHktYnktYXNzZW1ibHkuIEJlYXIgaW4gbWluZCB3ZSBoYXZlIHRvdGFsIGByIG5yb3codGhlX2RhdGEgJT4lIHNlbGVjdChhc3NlbWJseV9pZCkgJT4lIGRpc3RpbmN0KCkpDQpgIGFzc2VtYmxpZXMsIG9mIHdoaWNoIGByIG5yb3codGhlX2RhdGEgJT4lIGZpbHRlcihjb21tdW5pdHkgPT0gIk1HNWEiKSAlPiUgc2VsZWN0KGFzc2VtYmx5X2lkKSAlPiUgZGlzdGluY3QoKSkNCmAgYXJlIE1HNWEuDQoNCkxvb2tzIGxpa2UgdGhpcyB3aWxsIG5lZWQgc3BlY2lhbCB0cmVhdG1lbnQgd2l0aCBhIFNoaW55IGFwcCwgYmV5b25kIHRoZSBzY29wZSBvZiB0aGlzIGV4cGxvcmF0b3J5IG5vdGVib29rLg0KDQoNCg0KDQoNCg==